APIM - Control - Pushing Separate

knitr::opts_chunk$set(echo = TRUE)

cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)

system("shutdown /a")

## [1] 1116

# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
report_hurdle = TRUE
do_priorsense = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 2000 # 10'000 per chain to achieve 40'000
warmup = 1000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = as.character(iterations)

df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = FALSE, use_mi = use_mi)[[1]]

Constructing scales reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.3468 0.0000 5.0000 2323

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')

# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')

rows_to_pack <- list(
  "Within-Person Effects" = c(2,9),
  "Between-Person Effects" = c(10,16),
  "Random Effects" = c(17, 23), 
  "Additional Parameters" = c(24,24)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,9+5),
  "Between-Person Effects" = c(10+5,16+5),
  "Random Effects" = c(17+5, 23+5), 
  "Additional Parameters" = c(24+5,24+6)
  )

HURDLE MODELS

# For indistinguishable Dyads
model_rows_fixed_hu <- c(
    'Intercept', 
    'hu_Intercept',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    
    # HURDLE MODEL
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'hu_persuasion_self_cw', 
    'hu_persuasion_partner_cw', 
    'hu_pressure_self_cw', 
    'hu_pressure_partner_cw', 
    'hu_pushing_self_cw', 
    'hu_pushing_partner_cw', 
    'hu_day', 
    'hu_weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'hu_persuasion_self_cb',
    'hu_persuasion_partner_cb',
    'hu_pressure_self_cb',
    'hu_pressure_partner_cb',
    'hu_pushing_self_cb',
    'hu_pushing_partner_cb',
    'hu_weartime_self_cb'
  )

model_rows_fixed_hu_ordinal <- c(
  model_rows_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed_hu[3:length(model_rows_fixed_hu)]
)


model_rows_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(hu_Intercept)',
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # HURDLE
  'sd(hu_persuasion_self_cw)',
  'sd(hu_persuasion_partner_cw)',
  'sd(hu_pressure_self_cw)',
  'sd(hu_pressure_partner_cw)',
  'sd(hu_pushing_self_cw)',
  'sd(hu_pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_hu_ordinal <- c(model_rows_random_hu,'disc')

# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
    "Intercept", 
    "Hurdle Intercept",
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime",
    
    # HURDLE
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Hu Daily persuasion experienced", 
    "Hu Daily persuasion utilized (partner's view)", # OR partner received
    "Hu Daily pressure experienced", 
    "Hu Daily pressure utilized (partner's view)", 
    "Hu Daily pushing experienced", 
    "Hu Daily pushing utilized (partner's view)", 
    "Hu Day", 
    "Hu Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Hu Mean persuasion experienced", 
    "Hu Mean persuasion utilized (partner's view)", 
    "Hu Mean pressure experienced", 
    "Hu Mean pressure utilized (partner's view)", 
    "Hu Mean pushing experienced", 
    "Hu Mean pushing utilized (partner's view)", 
    "Hu Mean weartime"
  )



model_rownames_fixed_hu_ordinal <- c(
  model_rownames_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed_hu[3:length(model_rownames_fixed_hu)]
)



model_rownames_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Hurdle Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  
  # Hurdle
  "sd(Hu Daily persuasion experienced)", 
  "sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Hu Daily pressure experienced)", 
  "sd(Hu Daily pressure utilized (partner's view))", 
  "sd(Hu Daily pushing experienced)", 
  "sd(Hu Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_hu_ordinal <- c(model_rownames_random_hu,'disc')

rows_to_pack_hu <- list(
  "Conditional Within-Person Effects" = c(3,10),
  "Conditional Between-Person Effects" = c(11,17),
  
  "Hurdle Within-Person Effects" = c(18,25),
  "Hurdle Between-Person Effects" = c(26,32),
  
  "Random Effects" = c(33, 46), 
  "Additional Parameters" = c(47,47)
  )

rows_to_pack_hu_ordinal <- list(
  "Intercepts" = c(1,7),
  "Conditional Within-Person Effects" = c(3+5,10+5),
  "Conditional Between-Person Effects" = c(11+5,17+5),
  
  "Hurdle Within-Person Effects" = c(18+5,25+5),
  "Hurdle Between-Person Effects" = c(26+5,32+5),
  
  "Random Effects" = c(33+5, 46+5), 
  "Additional Parameters" = c(47+5,47+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 

## [1]   0 720

hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

WITHOUT PUSHING

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_Nopushing", suffix))
  #, file_refit = 'always'
)

## Warning: Rows containing NAs were excluded from the model.

pa_sub_digest <- digest::digest(pa_sub)

check_brms(pa_sub, log_pp_check = TRUE)

## 
## Divergences:

## 0 of 4000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 4000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -10536.2 166.2
## p_loo       144.3   4.8
## looic     21072.3 332.5
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.5]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3735  100.0%  512     
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.

DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 8, observations = 3736, p-value = 0.46
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000267666 0.003211991
## sample estimates:
## outlier frequency (expected: 0.00152837259100642 ) 
##                                        0.002141328

# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
  pa_sub, 
  stats_to_report = c('pd','ROPE'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 

## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.

# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub_continuous <- summarize_brms(
  pa_sub, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_continuous,
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 

## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.

# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous

columns_to_replace <- c("ROPE", "inside ROPE")

summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <- 
  summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack_hu)

	exp(Est.)	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	48.72***	3.08	[42.92, 55.10]	1.000	[0.92, 1.08]	0.000	1.002	964	1749
Hurdle Intercept	0.82	0.14	[ 0.59, 1.13]	0.882	[0.84, 1.20]	0.433	1.005	552	1254
Conditional Within-Person Effects
Daily persuasion experienced	1.03	0.03	[ 0.97, 1.09]	0.853	[0.92, 1.08]	0.965	1.001	1471	2332
Daily persuasion utilized (partner’s view)	1.03	0.02	[ 0.98, 1.08]	0.877	[0.92, 1.08]	0.987	1.002	2066	2263
Daily pressure experienced	0.90*	0.04	[ 0.81, 0.99]	0.983	[0.92, 1.08]	0.260	1.000	3707	2336
Daily pressure utilized (partner’s view)	0.94	0.04	[ 0.86, 1.02]	0.928	[0.92, 1.08]	0.634	1.002	4371	3030
Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Day	0.98	0.06	[ 0.87, 1.11]	0.623	[0.92, 1.08]	0.766	1.000	5745	2727
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Conditional Between-Person Effects
Mean persuasion experienced	1.11	0.15	[ 0.85, 1.47]	0.775	[0.92, 1.08]	0.336	1.004	817	1463
Mean persuasion utilized (partner’s view)	1.09	0.15	[ 0.83, 1.45]	0.737	[0.92, 1.08]	0.354	1.007	812	1214
Mean pressure experienced	1.12	0.18	[ 0.81, 1.55]	0.752	[0.92, 1.08]	0.287	1.005	1015	1888
Mean pressure utilized (partner’s view)	0.91	0.16	[ 0.64, 1.28]	0.708	[0.92, 1.08]	0.305	1.003	1145	1670
Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced	1.61***	0.10	[ 1.43, 1.84]	1.000	[0.84, 1.20]	0.000	1.001	1969	2771
Hu Daily persuasion utilized (partner’s view)	1.43***	0.08	[ 1.28, 1.61]	1.000	[0.84, 1.20]	0.000	1.001	3179	3039
Hu Daily pressure experienced	1.01	0.15	[ 0.76, 1.39]	0.525	[0.84, 1.20]	0.779	1.000	3515	2678
Hu Daily pressure utilized (partner’s view)	1.74***	0.32	[ 1.26, 2.83]	1.000	[0.84, 1.20]	0.009	1.001	2286	1405
Hu Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Day	0.92	0.12	[ 0.71, 1.20]	0.730	[0.84, 1.20]	0.759	1.001	7057	2959
Hu Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced	1.56	0.55	[ 0.76, 3.18]	0.889	[0.84, 1.20]	0.183	1.011	579	1175
Hu Mean persuasion utilized (partner’s view)	1.55	0.55	[ 0.76, 3.19]	0.881	[0.84, 1.20]	0.188	1.009	606	1142
Hu Mean pressure experienced	0.33**	0.14	[ 0.15, 0.77]	0.996	[0.84, 1.20]	0.013	1.007	782	1578
Hu Mean pressure utilized (partner’s view)	0.58	0.24	[ 0.25, 1.33]	0.911	[0.84, 1.20]	0.139	1.007	803	1296
Hu Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.32	0.04	[0.24, 0.42]	NA	NA	NA	1.008	985	1765
sd(Hurdle Intercept)	0.89	0.12	[0.70, 1.18]	NA	NA	NA	1.001	902	1729
sd(Daily persuasion experienced)	0.12	0.02	[0.08, 0.17]	NA	NA	NA	1.003	1813	2949
sd(Daily persuasion utilized (partner’s view))	0.09	0.02	[0.05, 0.14]	NA	NA	NA	1.006	2186	2589
sd(Daily pressure experienced)	0.08	0.07	[0.00, 0.25]	NA	NA	NA	1.001	1377	1876
sd(Daily pressure utilized (partner’s view))	0.07	0.05	[0.00, 0.19]	NA	NA	NA	1.003	1648	1853
sd(Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily persuasion experienced)	0.21	0.07	[0.08, 0.37]	NA	NA	NA	1.002	1341	977
sd(Hu Daily persuasion utilized (partner’s view))	0.19	0.07	[0.05, 0.34]	NA	NA	NA	1.001	1429	1309
sd(Hu Daily pressure experienced)	0.23	0.20	[0.01, 0.80]	NA	NA	NA	1.003	1209	1877
sd(Hu Daily pressure utilized (partner’s view))	0.31	0.25	[0.02, 1.01]	NA	NA	NA	1.004	1222	1709
sd(Hu Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	0.69	0.01	[0.67, 0.71]	NA	NA	NA	1.001	5060	3213

# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()

## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.77), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.79), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.77), b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.82), b_hu_pressure_self_cb and
##   b_hu_persuasion_self_cb (r = 0.81), b_hu_pressure_partner_cb and
##   b_hu_persuasion_partner_cb (r = 0.81). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)

## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.

print(conds_eff)

$persuasion_self_cw

## Warning: Removed 157 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 173 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Picking joint bandwidth of 0.011

$persuasion_partner_cw

## Warning: Removed 143 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 145 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Picking joint bandwidth of 0.00867

$pressure_self_cw

## Warning: Removed 46 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 55 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Picking joint bandwidth of 0.016

$pressure_partner_cw

## Warning: Removed 49 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 85 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Picking joint bandwidth of 0.0275

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 

## [1]   5.75 971.25

hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_nopushing", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

pa_obj_log_digest <- digest::digest(pa_obj_log)

check_brms(pa_obj_log, log_pp_check = TRUE)

## 
## Divergences:

## 0 of 4000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 4000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 4000 by 3337 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -18783.3  68.4
## p_loo        74.9   3.4
## looic     37566.6 136.7
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.7]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')

## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 27, observations = 3337, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.001498352 0.005394067
## sample estimates:
## outlier frequency (expected: 0.00329937069223854 ) 
##                                          0.0080911

# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)

## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.

	exp(Est.)	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	116.98***	6.11	[105.53, 130.15]	1.000	[0.94, 1.07]	0.000	1.003	583	1345
Within-Person Effects
Daily persuasion experienced	1.04*	0.01	[ 1.01, 1.07]	0.989	[0.94, 1.07]	0.975	1.001	2606	2982
Daily persuasion utilized (partner’s view)	1.02	0.02	[ 0.99, 1.06]	0.936	[0.94, 1.07]	0.995	1.000	3058	3179
Daily pressure experienced	0.95	0.03	[ 0.89, 1.02]	0.935	[0.94, 1.07]	0.678	1.000	4455	3085
Daily pressure utilized (partner’s view)	0.99	0.03	[ 0.93, 1.06]	0.627	[0.94, 1.07]	0.936	1.001	4236	2961
Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Day	0.97	0.03	[ 0.91, 1.04]	0.774	[0.94, 1.07]	0.856	1.001	7048	2972
Daily weartime	1.00***	0.00	[ 1.00, 1.00]	1.000	[0.94, 1.07]	1.000	1.001	3742	2490
Between-Person Effects
Mean persuasion experienced	1.11	0.14	[ 0.85, 1.42]	0.788	[0.94, 1.07]	0.279	1.007	529	1039
Mean persuasion utilized (partner’s view)	1.04	0.13	[ 0.79, 1.34]	0.611	[0.94, 1.07]	0.374	1.008	538	1075
Mean pressure experienced	0.89	0.13	[ 0.68, 1.20]	0.782	[0.94, 1.07]	0.244	1.004	669	1436
Mean pressure utilized (partner’s view)	1.02	0.14	[ 0.78, 1.37]	0.578	[0.94, 1.07]	0.362	1.005	643	1395
Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean weartime	1.00	0.00	[ 1.00, 1.00]	0.912	[0.94, 1.07]	1.000	1.000	4781	3605
Random Effects
sd(Intercept)	0.30	0.04	[0.23, 0.39]	NA	NA	NA	1.005	957	1591
sd(Daily persuasion experienced)	0.05	0.01	[0.03, 0.08]	NA	NA	NA	1.002	2335	2032
sd(Daily persuasion utilized (partner’s view))	0.06	0.02	[0.03, 0.09]	NA	NA	NA	1.001	1909	2207
sd(Daily pressure experienced)	0.05	0.04	[0.00, 0.14]	NA	NA	NA	1.001	1627	2057
sd(Daily pressure utilized (partner’s view))	0.04	0.03	[0.00, 0.12]	NA	NA	NA	1.000	2070	1630
sd(Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	0.58	0.01	[0.56, 0.59]	NA	NA	NA	1.004	7025	2589

plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()

## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.88), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.85), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.79), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.77), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.87). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)

## [1] 0 5

hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_nopushing", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

mood_gauss_digest <- digest::digest(mood_gauss)

check_brms(mood_gauss, log_pp_check = FALSE)

## 
## Divergences:

## 0 of 4000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 4000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -5189.1  59.2
## p_loo        68.2   3.0
## looic     10378.3 118.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3735  100.0%  750     
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.

DHARMa.check_brms.all(mood_gauss, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 32, observations = 3736, p-value =
## 2.24e-11
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.005865832 0.012070325
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                           0.00856531

summarize_brms(
  mood_gauss, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)

	Est.	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	3.70***	0.11	[ 3.49, 3.90]	1.000	[-0.11, 0.11]	0.000	1.006	417	729
Within-Person Effects
Daily persuasion experienced	0.01	0.02	[-0.03, 0.05]	0.663	[-0.11, 0.11]	1.000	1.001	3672	2684
Daily persuasion utilized (partner’s view)	0.04	0.02	[-0.01, 0.08]	0.938	[-0.11, 0.11]	0.999	1.000	3341	3254
Daily pressure experienced	-0.03	0.05	[-0.14, 0.07]	0.743	[-0.11, 0.11]	0.940	1.000	4889	2934
Daily pressure utilized (partner’s view)	0.01	0.05	[-0.10, 0.11]	0.559	[-0.11, 0.11]	0.965	1.001	3761	2704
Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Day	0.25***	0.06	[ 0.15, 0.36]	1.000	[-0.11, 0.11]	0.006	1.001	7106	2660
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean persuasion experienced	0.42	0.25	[-0.10, 0.91]	0.942	[-0.11, 0.11]	0.099	1.010	499	775
Mean persuasion utilized (partner’s view)	0.33	0.25	[-0.18, 0.82]	0.902	[-0.11, 0.11]	0.152	1.009	487	558
Mean pressure experienced	-0.36	0.26	[-0.91, 0.19]	0.903	[-0.11, 0.11]	0.139	1.007	571	954
Mean pressure utilized (partner’s view)	-0.29	0.27	[-0.82, 0.27]	0.851	[-0.11, 0.11]	0.180	1.007	549	835
Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.61	0.07	[0.48, 0.79]	NA	NA	NA	1.009	812	1382
sd(Daily persuasion experienced)	0.05	0.03	[0.00, 0.11]	NA	NA	NA	1.005	1000	1287
sd(Daily persuasion utilized (partner’s view))	0.08	0.03	[0.03, 0.14]	NA	NA	NA	1.003	1363	1122
sd(Daily pressure experienced)	0.06	0.06	[0.00, 0.23]	NA	NA	NA	1.001	1944	2092
sd(Daily pressure utilized (partner’s view))	0.07	0.06	[0.00, 0.25]	NA	NA	NA	1.002	2091	2255
sd(Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	0.96	0.01	[0.94, 0.98]	NA	NA	NA	1.000	8235	2720

plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()

## Possible multicollinearity between b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.86), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.86), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.82). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

Reactance

range(df_double$reactance, na.rm = T) 

## [1] 0 5

hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal_nopushing", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

reactance_ordinal_digest <- digest::digest(reactance_ordinal)

check_brms(reactance_ordinal)

## 
## Divergences:

## 0 of 4000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 4000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 5 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -683.2 31.6
## p_loo        59.6  4.5
## looic      1366.3 63.3
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     751   99.3%   224     
##    (0.7, 1]   (bad)        5    0.7%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 0.36
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.002645503
## sample estimates:
## outlier frequency (expected: 0.000304232804232804 ) 
##                                         0.001322751

summarize_brms(
  reactance_ordinal, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)

	OR	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercepts
Intercept	NA	NA	NA	NA	NA	NA	NA	NA	NA
Intercept[1]	3.48***	0.86	[ 2.16, 5.75]	1.000	[0.84, 1.20]	0.000	1.000	3268	2969
Intercept[2]	7.32***	1.90	[ 4.48, 12.47]	1.000	[0.84, 1.20]	0.000	1.001	3373	3061
Intercept[3]	19.78***	5.57	[ 11.49, 35.14]	1.000	[0.84, 1.20]	0.000	1.001	3544	3388
Intercept[4]	84.24***	29.47	[ 45.16, 173.12]	1.000	[0.84, 1.20]	0.000	1.000	4264	3330
Intercept[5]	2760.30***	1715.32	[929.55, 10049.83]	1.000	[0.84, 1.20]	0.000	1.000	5669	2968
Within-Person Effects
Daily persuasion experienced	0.85*	0.07	[ 0.72, 0.99]	0.979	[0.84, 1.20]	0.596	1.000	4158	3111
Daily persuasion utilized (partner’s view)	1.01	0.09	[ 0.84, 1.20]	0.548	[0.84, 1.20]	0.950	1.001	4250	2828
Daily pressure experienced	1.98**	0.34	[ 1.31, 2.76]	0.999	[0.84, 1.20]	0.011	1.001	2901	3115
Daily pressure utilized (partner’s view)	1.15	0.23	[ 0.72, 1.85]	0.749	[0.84, 1.20]	0.509	1.001	3088	2320
Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Day	1.63	0.55	[ 0.87, 3.10]	0.929	[0.84, 1.20]	0.161	1.000	6064	2790
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean persuasion experienced	1.08	0.50	[ 0.42, 2.76]	0.567	[0.84, 1.20]	0.309	1.005	1276	2066
Mean persuasion utilized (partner’s view)	0.79	0.42	[ 0.29, 2.26]	0.665	[0.84, 1.20]	0.243	1.004	1346	2146
Mean pressure experienced	3.30*	1.70	[ 1.21, 9.23]	0.990	[0.84, 1.20]	0.019	1.004	1610	2466
Mean pressure utilized (partner’s view)	1.23	0.70	[ 0.41, 3.54]	0.637	[0.84, 1.20]	0.226	1.003	1666	2712
Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.78	0.19	[0.46, 1.22]	NA	NA	NA	1.001	1742	2619
sd(Daily persuasion experienced)	0.16	0.12	[0.01, 0.44]	NA	NA	NA	1.008	749	1628
sd(Daily persuasion utilized (partner’s view))	0.18	0.13	[0.01, 0.45]	NA	NA	NA	1.002	1314	1689
sd(Daily pressure experienced)	0.47	0.22	[0.07, 1.02]	NA	NA	NA	1.007	1083	919
sd(Daily pressure utilized (partner’s view))	0.31	0.31	[0.02, 1.37]	NA	NA	NA	1.002	1028	1907
sd(Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	NA	NA	NA	NA	NA	NA	NA	NA	NA
disc	1.00	0.00	[1.00, 1.00]	NA	NA	NA	NA	NA	NA

plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()

## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.75), b_Intercept[4] and b_Intercept[3] (r = 0.83),
##   b_pressure_self_cb and b_persuasion_self_cb (r = 0.79),
##   b_pressure_partner_cb and b_persuasion_partner_cb (r = 0.84). This might
##   lead to inappropriate results. See 'Details' in '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="03_SensitivityPushingSeparate_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_nopushing_", suffix))
  #, file_refit = 'always'
)

## Warning: Rows containing NAs were excluded from the model.

is_reactance_digest <- digest::digest(is_reactance)

check_brms(is_reactance)

## 
## Divergences:

## 0 of 4000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 4000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 31 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -361.7 15.6
## p_loo        63.7  4.9
## looic       723.3 31.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.2, 1.5]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     725   95.9%   78      
##    (0.7, 1]   (bad)       30    4.0%   <NA>    
##    (1, Inf)   (very bad)   1    0.1%   <NA>    
## See help('pareto-k-diagnostic') for details.

DHARMa.check_brms.all(is_reactance, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 2, observations = 756, p-value = 0.6663
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0003205434 0.0095234999
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.002645503

summarize_brms(
  is_reactance, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)

	OR	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	0.32***	0.08	[0.18, 0.53]	1.000	[0.83, 1.20]	0.000	1.002	2308	2705
Within-Person Effects
Daily persuasion experienced	0.86	0.08	[0.71, 1.01]	0.964	[0.83, 1.20]	0.612	1.000	2887	2614
Daily persuasion utilized (partner’s view)	1.10	0.14	[0.86, 1.45]	0.777	[0.83, 1.20]	0.729	1.001	2245	2471
Daily pressure experienced	2.09*	0.58	[1.19, 4.20]	0.991	[0.83, 1.20]	0.024	1.001	1589	1427
Daily pressure utilized (partner’s view)	1.34	0.53	[0.57, 3.87]	0.780	[0.83, 1.20]	0.276	1.001	1658	1539
Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Day	1.86	0.67	[0.89, 3.88]	0.955	[0.83, 1.20]	0.092	1.001	4640	2592
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean persuasion experienced	1.86	1.07	[0.58, 5.93]	0.859	[0.83, 1.20]	0.144	1.005	1143	2107
Mean persuasion utilized (partner’s view)	1.03	0.66	[0.31, 3.66]	0.520	[0.83, 1.20]	0.233	1.004	1265	2190
Mean pressure experienced	17.76***	16.40	[3.37, 123.01]	1.000	[0.83, 1.20]	0.000	1.002	1423	2478
Mean pressure utilized (partner’s view)	1.21	1.21	[0.17, 7.77]	0.574	[0.83, 1.20]	0.136	1.001	1712	2485
Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	1.21	0.25	[0.80, 1.78]	NA	NA	NA	1.001	1396	1835
sd(Daily persuasion experienced)	0.19	0.13	[0.01, 0.48]	NA	NA	NA	1.004	958	1649
sd(Daily persuasion utilized (partner’s view))	0.43	0.17	[0.12, 0.83]	NA	NA	NA	1.003	1279	1064
sd(Daily pressure experienced)	0.79	0.53	[0.07, 2.12]	NA	NA	NA	1.004	481	715
sd(Daily pressure utilized (partner’s view))	0.74	0.64	[0.04, 2.71]	NA	NA	NA	1.002	1030	1680
sd(Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	NA	NA	NA	NA	NA	NA	NA	NA	NA

plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

Report All Models

process_model_component <- function(obj) {
  # Convert to string, modify, and evaluate
  name <- deparse(substitute(obj))
  if (report_hurdle) name <- paste0(name, '_hu')
  if (report_ordinal) name <- paste0(name, '_ordinal')
  return(get(name, envir = parent.frame()))
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI', 'pd'),
  
  model_rows_random = process_model_component(model_rows_random),
  model_rows_fixed = process_model_component(model_rows_fixed),
  model_rownames_random = process_model_component(model_rownames_random),
  model_rownames_fixed = process_model_component(model_rownames_fixed)
) 

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”

# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),  
      "Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5), 
      "Mood Gaussian" = (length(all_models) / 5),
      "Reactance Ordinal" = (length(all_models) / 5),
      "Reactance Dichotome" = (length(all_models) / 5))
  )


export_xlsx(
  summary_all_models, 
  rows_to_pack = process_model_component(rows_to_pack),
  file.path("Output", "AllModelsNoPushing.xlsx"), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
print(summary_all_models)

	Subjective MVPA Hurdle Lognormal			Device-Based MVPA Log (Gaussian)			Mood Gaussian			Reactance Ordinal			Reactance Dichotome
	exp(Est.) pa_sub	95% CI pa_sub	pd pa_sub	exp(Est.) pa_obj_log	95% CI pa_obj_log	pd pa_obj_log	Est. mood_gauss	95% CI mood_gauss	pd mood_gauss	OR reactance_ordinal	95% CI reactance_ordinal	pd reactance_ordinal	OR is_reactance	95% CI is_reactance	pd is_reactance
Intercept	48.72***	[42.92, 55.10]	1.000	116.98***	[105.53, 130.15]	1.000	3.70***	[ 3.49, 3.90]	1.000	NA	NA	NA	0.32***	[0.18, 0.53]	1.000
Hurdle Intercept	0.82	[ 0.59, 1.13]	0.882	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Conditional Within-Person Effects
Daily persuasion experienced	1.03	[ 0.97, 1.09]	0.853	1.04*	[ 1.01, 1.07]	0.989	0.01	[-0.03, 0.05]	0.663	0.85*	[ 0.72, 0.99]	0.979	0.86	[0.71, 1.01]	0.964
Daily persuasion utilized (partner’s view)	1.03	[ 0.98, 1.08]	0.877	1.02	[ 0.99, 1.06]	0.936	0.04	[-0.01, 0.08]	0.938	1.01	[ 0.84, 1.20]	0.548	1.10	[0.86, 1.45]	0.777
Daily pressure experienced	0.90*	[ 0.81, 0.99]	0.983	0.95	[ 0.89, 1.02]	0.935	-0.03	[-0.14, 0.07]	0.743	1.98**	[ 1.31, 2.76]	0.999	2.09*	[1.19, 4.20]	0.991
Daily pressure utilized (partner’s view)	0.94	[ 0.86, 1.02]	0.928	0.99	[ 0.93, 1.06]	0.627	0.01	[-0.10, 0.11]	0.559	1.15	[ 0.72, 1.85]	0.749	1.34	[0.57, 3.87]	0.780
Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Day	0.98	[ 0.87, 1.11]	0.623	0.97	[ 0.91, 1.04]	0.774	0.25***	[ 0.15, 0.36]	1.000	1.63	[ 0.87, 3.10]	0.929	1.86	[0.89, 3.88]	0.955
Daily weartime	NA	NA	NA	1.00***	[ 1.00, 1.00]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA
Conditional Between-Person Effects
Mean persuasion experienced	1.11	[ 0.85, 1.47]	0.775	1.11	[ 0.85, 1.42]	0.788	0.42	[-0.10, 0.91]	0.942	1.08	[ 0.42, 2.76]	0.567	1.86	[0.58, 5.93]	0.859
Mean persuasion utilized (partner’s view)	1.09	[ 0.83, 1.45]	0.737	1.04	[ 0.79, 1.34]	0.611	0.33	[-0.18, 0.82]	0.902	0.79	[ 0.29, 2.26]	0.665	1.03	[0.31, 3.66]	0.520
Mean pressure experienced	1.12	[ 0.81, 1.55]	0.752	0.89	[ 0.68, 1.20]	0.782	-0.36	[-0.91, 0.19]	0.903	3.30*	[ 1.21, 9.23]	0.990	17.76***	[3.37, 123.01]	1.000
Mean pressure utilized (partner’s view)	0.91	[ 0.64, 1.28]	0.708	1.02	[ 0.78, 1.37]	0.578	-0.29	[-0.82, 0.27]	0.851	1.23	[ 0.41, 3.54]	0.637	1.21	[0.17, 7.77]	0.574
Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean weartime	NA	NA	NA	1.00	[ 1.00, 1.00]	0.912	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced	1.61***	[ 1.43, 1.84]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily persuasion utilized (partner’s view)	1.43***	[ 1.28, 1.61]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pressure experienced	1.01	[ 0.76, 1.39]	0.525	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pressure utilized (partner’s view)	1.74***	[ 1.26, 2.83]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Day	0.92	[ 0.71, 1.20]	0.730	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced	1.56	[ 0.76, 3.18]	0.889	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean persuasion utilized (partner’s view)	1.55	[ 0.76, 3.19]	0.881	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pressure experienced	0.33**	[ 0.15, 0.77]	0.996	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pressure utilized (partner’s view)	0.58	[ 0.25, 1.33]	0.911	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.32	[0.24, 0.42]	NA	0.30	[0.23, 0.39]	NA	0.61	[0.48, 0.79]	NA	0.78	[0.46, 1.22]	NA	1.21	[0.80, 1.78]	NA
sd(Hurdle Intercept)	0.89	[0.70, 1.18]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion experienced)	0.12	[0.08, 0.17]	NA	0.05	[0.03, 0.08]	NA	0.05	[0.00, 0.11]	NA	0.16	[0.01, 0.44]	NA	0.19	[0.01, 0.48]	NA
sd(Daily persuasion utilized (partner’s view))	0.09	[0.05, 0.14]	NA	0.06	[0.03, 0.09]	NA	0.08	[0.03, 0.14]	NA	0.18	[0.01, 0.45]	NA	0.43	[0.12, 0.83]	NA
sd(Daily pressure experienced)	0.08	[0.00, 0.25]	NA	0.05	[0.00, 0.14]	NA	0.06	[0.00, 0.23]	NA	0.47	[0.07, 1.02]	NA	0.79	[0.07, 2.12]	NA
sd(Daily pressure utilized (partner’s view))	0.07	[0.00, 0.19]	NA	0.04	[0.00, 0.12]	NA	0.07	[0.00, 0.25]	NA	0.31	[0.02, 1.37]	NA	0.74	[0.04, 2.71]	NA
sd(Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily persuasion experienced)	0.21	[0.08, 0.37]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily persuasion utilized (partner’s view))	0.19	[0.05, 0.34]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pressure experienced)	0.23	[0.01, 0.80]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pressure utilized (partner’s view))	0.31	[0.02, 1.01]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	0.69	[0.67, 0.71]	NA	0.58	[0.56, 0.59]	NA	0.96	[0.94, 0.98]	NA	NA	NA	NA	NA	NA	NA

ONLY PUSHING

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_onlypushing", suffix))
  #, file_refit = 'always'
)

## Warning: Rows containing NAs were excluded from the model.

## Start sampling

pa_sub_digest <- digest::digest(pa_sub)

check_brms(pa_sub, log_pp_check = TRUE)

## 
## Divergences:

## 0 of 4000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 4000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
## Computed from 4000 by 1342 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -5692.9  78.2
## p_loo        97.8   4.1
## looic     11385.8 156.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1341  99.9%   633     
##    (0.7, 1]   (bad)         1   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 4, observations = 1342, p-value = 0.26
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.003371833
## sample estimates:
## outlier frequency (expected: 0.00149776453055142 ) 
##                                        0.002980626

# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
  pa_sub, 
  stats_to_report = c('pd','ROPE'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 

## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.

# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub_continuous <- summarize_brms(
  pa_sub, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_continuous,
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 

## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.

# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous

columns_to_replace <- c("ROPE", "inside ROPE")

summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <- 
  summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack_hu)

	exp(Est.)	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	57.28***	4.53	[49.10, 66.97]	1.000	[0.93, 1.08]	0.000	1.002	1327	2040
Hurdle Intercept	3.79***	0.98	[ 2.27, 6.26]	1.000	[0.84, 1.20]	0.000	1.005	1088	2463
Conditional Within-Person Effects
Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing experienced	0.97	0.03	[ 0.91, 1.03]	0.856	[0.93, 1.08]	0.931	1.002	4427	2662
Daily pushing utilized (partner’s view)	0.97	0.03	[ 0.90, 1.03]	0.847	[0.93, 1.08]	0.907	1.001	3679	2873
Day	0.91	0.07	[ 0.79, 1.06]	0.894	[0.93, 1.08]	0.398	1.000	8744	2600
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Conditional Between-Person Effects
Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing experienced	1.32*	0.18	[ 1.00, 1.71]	0.977	[0.93, 1.08]	0.060	1.005	1217	2176
Mean pushing utilized (partner’s view)	1.18	0.15	[ 0.90, 1.53]	0.883	[0.93, 1.08]	0.222	1.005	1265	2076
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pushing experienced	1.55**	0.23	[ 1.17, 2.24]	0.998	[0.84, 1.20]	0.035	1.002	3127	2531
Hu Daily pushing utilized (partner’s view)	1.84***	0.30	[ 1.38, 2.68]	1.000	[0.84, 1.20]	0.002	1.001	2766	2989
Hu Day	0.61	0.15	[ 0.37, 1.00]	0.974	[0.84, 1.20]	0.100	1.003	7632	2768
Hu Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pushing experienced	0.46**	0.14	[ 0.26, 0.81]	0.997	[0.84, 1.20]	0.018	1.001	2675	2824
Hu Mean pushing utilized (partner’s view)	0.62	0.18	[ 0.33, 1.08]	0.952	[0.84, 1.20]	0.136	1.001	2568	2699
Hu Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.37	0.06	[0.29, 0.50]	NA	NA	NA	1.002	1283	2370
sd(Hurdle Intercept)	1.25	0.19	[0.93, 1.74]	NA	NA	NA	1.002	1448	2426
sd(Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing experienced)	0.06	0.04	[0.00, 0.16]	NA	NA	NA	1.001	870	1486
sd(Daily pushing utilized (partner’s view))	0.09	0.04	[0.01, 0.18]	NA	NA	NA	1.004	977	1142
sd(Hu Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pushing experienced)	0.46	0.24	[0.07, 1.02]	NA	NA	NA	1.000	947	1689
sd(Hu Daily pushing utilized (partner’s view))	0.51	0.24	[0.10, 1.13]	NA	NA	NA	1.001	1088	1582
Additional Parameters
sigma	0.67	0.02	[0.64, 0.71]	NA	NA	NA	1.002	6374	2773

# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()

## Possible multicollinearity between b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.77). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)

## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.

print(conds_eff)

$pushing_self_cw

## Warning: Removed 200 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 200 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Picking joint bandwidth of 0.0172

$pushing_partner_cw

## Warning: Removed 130 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 180 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Picking joint bandwidth of 0.022

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 

## [1]   5.75 971.25

hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
   
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_onlypushing", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

pa_obj_log_digest <- digest::digest(pa_obj_log)

check_brms(pa_obj_log, log_pp_check = TRUE)

## 
## Divergences:

## 0 of 4000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 4000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 4000 by 1214 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -6914.2 39.4
## p_loo        58.6  4.3
## looic     13828.4 78.8
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')

## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 8, observations = 1214, p-value = 0.16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0008237232 0.0074135091
## sample estimates:
## outlier frequency (expected: 0.00339373970345964 ) 
##                                        0.006589786

# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)

## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.

	exp(Est.)	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	124.24***	7.52	[109.69, 139.58]	1.000	[0.94, 1.06]	0.000	1.002	967	1634
Within-Person Effects
Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing experienced	1.03	0.03	[ 0.97, 1.08]	0.807	[0.94, 1.06]	0.914	1.001	3221	2733
Daily pushing utilized (partner’s view)	1.02	0.02	[ 0.97, 1.06]	0.797	[0.94, 1.06]	0.979	1.000	5065	2748
Day	0.97	0.05	[ 0.87, 1.08]	0.718	[0.94, 1.06]	0.654	1.000	7350	3336
Daily weartime	1.00**	0.00	[ 1.00, 1.00]	0.998	[0.94, 1.06]	1.000	1.001	4680	2898
Between-Person Effects
Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing experienced	0.97	0.06	[ 0.86, 1.10]	0.680	[0.94, 1.06]	0.635	1.001	2223	2645
Mean pushing utilized (partner’s view)	1.04	0.07	[ 0.92, 1.18]	0.739	[0.94, 1.06]	0.578	1.002	2058	2692
Mean weartime	1.00	0.00	[ 1.00, 1.00]	0.812	[0.94, 1.06]	1.000	1.000	3879	3246
Random Effects
sd(Intercept)	0.31	0.04	[0.24, 0.42]	NA	NA	NA	1.005	890	1790
sd(Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing experienced)	0.09	0.04	[0.02, 0.17]	NA	NA	NA	1.004	887	821
sd(Daily pushing utilized (partner’s view))	0.03	0.03	[0.00, 0.10]	NA	NA	NA	1.001	1398	1731
Additional Parameters
sigma	0.54	0.01	[0.52, 0.57]	NA	NA	NA	1.001	6040	3037

plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)

## [1] 0 5

hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_onlypushing", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

mood_gauss_digest <- digest::digest(mood_gauss)

check_brms(mood_gauss, log_pp_check = FALSE)

## 
## Divergences:

## 0 of 4000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 4000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 4000 by 1342 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1798.8 39.9
## p_loo        50.3  3.0
## looic      3597.5 79.8
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.8]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

DHARMa.check_brms.all(mood_gauss, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 13, observations = 1342, p-value =
## 4.837e-06
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.005167734 0.016508158
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.009687034

summarize_brms(
  mood_gauss, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)

	Est.	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	3.78***	0.12	[ 3.56, 4.00]	1.000	[-0.11, 0.11]	0.000	1.013	560	1111
Within-Person Effects
Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing experienced	0.02	0.04	[-0.05, 0.08]	0.670	[-0.11, 0.11]	0.995	1.002	2823	2591
Daily pushing utilized (partner’s view)	0.09*	0.04	[ 0.01, 0.17]	0.984	[-0.11, 0.11]	0.715	1.002	1597	2207
Day	0.28***	0.09	[ 0.11, 0.46]	1.000	[-0.11, 0.11]	0.023	1.001	3871	3088
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing experienced	0.12	0.11	[-0.08, 0.32]	0.880	[-0.11, 0.11]	0.433	1.003	1269	1954
Mean pushing utilized (partner’s view)	0.08	0.11	[-0.12, 0.30]	0.786	[-0.11, 0.11]	0.554	1.002	1216	1958
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.61	0.08	[0.48, 0.79]	NA	NA	NA	1.011	648	1448
sd(Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing experienced)	0.07	0.04	[0.00, 0.16]	NA	NA	NA	1.000	1409	1126
sd(Daily pushing utilized (partner’s view))	0.14	0.04	[0.06, 0.24]	NA	NA	NA	1.002	1329	828
Additional Parameters
sigma	0.91	0.02	[0.87, 0.94]	NA	NA	NA	1.001	3887	2536

plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

range(df_double$reactance, na.rm = T) 

## [1] 0 5

hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal_onlypushing", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

reactance_ordinal_digest <- digest::digest(reactance_ordinal)

check_brms(reactance_ordinal)

## 
## Divergences:

## 0 of 4000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 4000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 4000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -363.6 23.5
## p_loo        33.6  2.9
## looic       727.1 47.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.6]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 0, observations = 403, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.00248139
## sample estimates:
## outlier frequency (expected: 0.000297766749379653 ) 
##                                                   0

summarize_brms(
  reactance_ordinal, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)

	OR	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercepts
Intercept	NA	NA	NA	NA	NA	NA	NA	NA	NA
Intercept[1]	5.44***	1.57	[ 3.06, 10.06]	1.000	[0.84, 1.20]	0.000	1.001	3521	2882
Intercept[2]	9.89***	3.14	[ 5.45, 18.74]	1.000	[0.84, 1.20]	0.000	1.002	3598	2951
Intercept[3]	23.10***	7.89	[ 12.38, 46.79]	1.000	[0.84, 1.20]	0.000	1.001	3933	2899
Intercept[4]	81.88***	35.09	[ 38.53, 199.63]	1.000	[0.84, 1.20]	0.000	1.000	4569	3054
Intercept[5]	771.25***	625.85	[202.52, 4855.18]	1.000	[0.84, 1.20]	0.000	1.001	5246	2545
Within-Person Effects
Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing experienced	1.31*	0.14	[ 1.04, 1.62]	0.987	[0.84, 1.20]	0.214	1.001	3641	2884
Daily pushing utilized (partner’s view)	0.96	0.13	[ 0.75, 1.25]	0.615	[0.84, 1.20]	0.810	1.002	4044	2794
Day	1.60	0.73	[ 0.68, 3.95]	0.854	[0.84, 1.20]	0.182	1.002	4822	3245
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing experienced	2.32**	0.65	[ 1.35, 3.99]	1.000	[0.84, 1.20]	0.005	1.001	3394	2987
Mean pushing utilized (partner’s view)	1.25	0.35	[ 0.70, 2.16]	0.783	[0.84, 1.20]	0.365	1.001	3738	2562
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.85	0.21	[0.49, 1.34]	NA	NA	NA	1.001	1569	2204
sd(Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing experienced)	0.17	0.14	[0.01, 0.50]	NA	NA	NA	1.004	1241	1651
sd(Daily pushing utilized (partner’s view))	0.16	0.14	[0.01, 0.57]	NA	NA	NA	1.003	1626	1714
Additional Parameters
sigma	NA	NA	NA	NA	NA	NA	NA	NA	NA
disc	1.00	0.00	[1.00, 1.00]	NA	NA	NA	NA	NA	NA

plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()

## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.75), b_Intercept[4] and b_Intercept[3] (r = 0.82). This might lead
##   to inappropriate results. See 'Details' in '?rope'.

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_onlypushing", suffix))
  #, file_refit = 'always'
)

## Warning: Rows containing NAs were excluded from the model.

is_reactance_digest <- digest::digest(is_reactance)

check_brms(is_reactance)

## 
## Divergences:

## 0 of 4000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 4000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 3 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
## Computed from 4000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -216.1 11.0
## p_loo        39.3  3.0
## looic       432.3 22.1
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     400   99.3%   458     
##    (0.7, 1]   (bad)        3    0.7%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

DHARMa.check_brms.all(is_reactance, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 403, p-value = 0.5534
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  6.282137e-05 1.374725e-02
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                           0.00248139

summarize_brms(
  is_reactance, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)

	OR	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	0.40**	0.11	[0.23, 0.68]	1.000	[0.83, 1.20]	0.002	1.000	3095	2830
Within-Person Effects
Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing experienced	1.37*	0.18	[1.06, 1.87]	0.993	[0.83, 1.20]	0.154	1.000	3238	2410
Daily pushing utilized (partner’s view)	1.16	0.26	[0.80, 2.12]	0.769	[0.83, 1.20]	0.509	1.000	2310	2170
Day	1.78	0.89	[0.69, 4.81]	0.874	[0.83, 1.20]	0.160	1.000	6984	2927
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing experienced	5.39***	2.58	[2.38, 14.94]	1.000	[0.83, 1.20]	0.000	1.001	2293	2628
Mean pushing utilized (partner’s view)	2.37	1.15	[0.97, 6.82]	0.972	[0.83, 1.20]	0.057	1.001	2556	2640
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	1.52	0.31	[1.03, 2.24]	NA	NA	NA	1.001	1431	2494
sd(Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing experienced)	0.18	0.15	[0.01, 0.56]	NA	NA	NA	1.003	1256	1876
sd(Daily pushing utilized (partner’s view))	0.42	0.27	[0.03, 1.15]	NA	NA	NA	1.002	1292	1879
Additional Parameters
sigma	NA	NA	NA	NA	NA	NA	NA	NA	NA

plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_self_cw

Report All Models

process_model_component <- function(obj) {
  # Convert to string, modify, and evaluate
  name <- deparse(substitute(obj))
  if (report_hurdle) name <- paste0(name, '_hu')
  if (report_ordinal) name <- paste0(name, '_ordinal')
  return(get(name, envir = parent.frame()))
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI', 'pd'),
  
  model_rows_random = process_model_component(model_rows_random),
  model_rows_fixed = process_model_component(model_rows_fixed),
  model_rownames_random = process_model_component(model_rownames_random),
  model_rownames_fixed = process_model_component(model_rownames_fixed)
) 

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”

# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),  
      "Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5), 
      "Mood Gaussian" = (length(all_models) / 5),
      "Reactance Ordinal" = (length(all_models) / 5),
      "Reactance Dichotome" = (length(all_models) / 5))
  )


export_xlsx(
  summary_all_models, 
  rows_to_pack = process_model_component(rows_to_pack),
  file.path("Output", "AllModelsOnlyPushing.xlsx"), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
print(summary_all_models)

	Subjective MVPA Hurdle Lognormal			Device-Based MVPA Log (Gaussian)			Mood Gaussian			Reactance Ordinal			Reactance Dichotome
	exp(Est.) pa_sub	95% CI pa_sub	pd pa_sub	exp(Est.) pa_obj_log	95% CI pa_obj_log	pd pa_obj_log	Est. mood_gauss	95% CI mood_gauss	pd mood_gauss	OR reactance_ordinal	95% CI reactance_ordinal	pd reactance_ordinal	OR is_reactance	95% CI is_reactance	pd is_reactance
Intercept	57.28***	[49.10, 66.97]	1.000	124.24***	[109.69, 139.58]	1.000	3.78***	[ 3.56, 4.00]	1.000	NA	NA	NA	0.40**	[0.23, 0.68]	1.000
Hurdle Intercept	3.79***	[ 2.27, 6.26]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Conditional Within-Person Effects
Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing experienced	0.97	[ 0.91, 1.03]	0.856	1.03	[ 0.97, 1.08]	0.807	0.02	[-0.05, 0.08]	0.670	1.31*	[ 1.04, 1.62]	0.987	1.37*	[1.06, 1.87]	0.993
Daily pushing utilized (partner’s view)	0.97	[ 0.90, 1.03]	0.847	1.02	[ 0.97, 1.06]	0.797	0.09*	[ 0.01, 0.17]	0.984	0.96	[ 0.75, 1.25]	0.615	1.16	[0.80, 2.12]	0.769
Day	0.91	[ 0.79, 1.06]	0.894	0.97	[ 0.87, 1.08]	0.718	0.28***	[ 0.11, 0.46]	1.000	1.60	[ 0.68, 3.95]	0.854	1.78	[0.69, 4.81]	0.874
Daily weartime	NA	NA	NA	1.00**	[ 1.00, 1.00]	0.998	NA	NA	NA	NA	NA	NA	NA	NA	NA
Conditional Between-Person Effects
Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing experienced	1.32*	[ 1.00, 1.71]	0.977	0.97	[ 0.86, 1.10]	0.680	0.12	[-0.08, 0.32]	0.880	2.32**	[ 1.35, 3.99]	1.000	5.39***	[2.38, 14.94]	1.000
Mean pushing utilized (partner’s view)	1.18	[ 0.90, 1.53]	0.883	1.04	[ 0.92, 1.18]	0.739	0.08	[-0.12, 0.30]	0.786	1.25	[ 0.70, 2.16]	0.783	2.37	[0.97, 6.82]	0.972
Mean weartime	NA	NA	NA	1.00	[ 1.00, 1.00]	0.812	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pushing experienced	1.55**	[ 1.17, 2.24]	0.998	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pushing utilized (partner’s view)	1.84***	[ 1.38, 2.68]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Day	0.61	[ 0.37, 1.00]	0.974	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pushing experienced	0.46**	[ 0.26, 0.81]	0.997	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pushing utilized (partner’s view)	0.62	[ 0.33, 1.08]	0.952	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.37	[0.29, 0.50]	NA	0.31	[0.24, 0.42]	NA	0.61	[0.48, 0.79]	NA	0.85	[0.49, 1.34]	NA	1.52	[1.03, 2.24]	NA
sd(Hurdle Intercept)	1.25	[0.93, 1.74]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing experienced)	0.06	[0.00, 0.16]	NA	0.09	[0.02, 0.17]	NA	0.07	[0.00, 0.16]	NA	0.17	[0.01, 0.50]	NA	0.18	[0.01, 0.56]	NA
sd(Daily pushing utilized (partner’s view))	0.09	[0.01, 0.18]	NA	0.03	[0.00, 0.10]	NA	0.14	[0.06, 0.24]	NA	0.16	[0.01, 0.57]	NA	0.42	[0.03, 1.15]	NA
sd(Hu Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pushing experienced)	0.46	[0.07, 1.02]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pushing utilized (partner’s view))	0.51	[0.10, 1.13]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	0.67	[0.64, 0.71]	NA	0.54	[0.52, 0.57]	NA	0.91	[0.87, 0.94]	NA	NA	NA	NA	NA	NA	NA

report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::report_packages()

• rmarkdown (version 2.29; Allaire J et al., 2024)
• beepr (version 2.0; Bååth R, 2024)
• R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
• R.oo (version 1.27.0; Bengtsson H, 2003)
• R.utils (version 2.12.3; Bengtsson H, 2023)
• brms (version 2.22.0; Bürkner P, 2017)
• posterior (version 1.6.0; Bürkner P et al., 2024)
• extrafont (version 0.19; Chang W, 2023)
• digest (version 0.6.37; Eddelbuettel D, 2024)
• Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
• bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
• lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
• DHARMa (version 0.4.7; Hartig F, 2024)
• fs (version 1.6.5; Hester J et al., 2024)
• wbCorr (version 0.1.22; Küng P, 2023)
• see (version 0.9.0; Lüdecke D et al., 2021)
• tibble (version 3.2.1; Müller K, Wickham H, 2023)
• patchwork (version 1.3.0; Pedersen T, 2024)
• R (version 4.4.1; R Core Team, 2024)
• openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
• ggplot2 (version 3.5.1; Wickham H, 2016)
• forcats (version 1.0.0; Wickham H, 2023)
• stringr (version 1.5.1; Wickham H, 2023)
• rvest (version 1.0.4; Wickham H, 2024)
• tidyverse (version 2.0.0; Wickham H et al., 2019)
• readxl (version 1.4.3; Wickham H, Bryan J, 2023)
• dplyr (version 1.1.4; Wickham H et al., 2023)
• purrr (version 1.0.2; Wickham H, Henry L, 2023)
• readr (version 2.1.5; Wickham H et al., 2024)
• xml2 (version 1.3.6; Wickham H et al., 2023)
• tidyr (version 1.3.1; Wickham H et al., 2024)
• ggridges (version 0.5.6; Wilke C, 2024)
• knitr (version 1.49; Xie Y, 2024)
• kableExtra (version 1.4.0; Zhu H, 2024)

report::cite_packages()

• Allaire J, Xie Y, Dervieux C, McPherson J, Luraschi J, Ushey K, Atkins A, Wickham H, Cheng J, Chang W, Iannone R (2024). rmarkdown: Dynamic Documents for R. R package version 2.29, https://github.com/rstudio/rmarkdown. Xie Y, Allaire J, Grolemund G (2018). R Markdown: The Definitive Guide. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9781138359338, https://bookdown.org/yihui/rmarkdown. Xie Y, Dervieux C, Riederer E (2020). R Markdown Cookbook. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9780367563837, https://bookdown.org/yihui/rmarkdown-cookbook.
• Bååth R (2024). beepr: Easily Play Notification Sounds on any Platform. R package version 2.0, https://CRAN.R-project.org/package=beepr.
• Bengtsson H (2003). “The R.oo package - Object-Oriented Programming with References Using Standard R Code.” In Hornik K, Leisch F, Zeileis A (eds.), Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003). https://www.r-project.org/conferences/DSC-2003/Proceedings/Bengtsson.pdf.
• Bengtsson H (2003). “The R.oo package - Object-Oriented Programming with References Using Standard R Code.” In Hornik K, Leisch F, Zeileis A (eds.), Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003). https://www.r-project.org/conferences/DSC-2003/Proceedings/Bengtsson.pdf.
• Bengtsson H (2023). R.utils: Various Programming Utilities. R package version 2.12.3, https://CRAN.R-project.org/package=R.utils.
• Bürkner P (2017). “brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software, 80(1), 1-28. doi:10.18637/jss.v080.i01 https://doi.org/10.18637/jss.v080.i01. Bürkner P (2018). “Advanced Bayesian Multilevel Modeling with the R Package brms.” The R Journal, 10(1), 395-411. doi:10.32614/RJ-2018-017 https://doi.org/10.32614/RJ-2018-017. Bürkner P (2021). “Bayesian Item Response Modeling in R with brms and Stan.” Journal of Statistical Software, 100(5), 1-54. doi:10.18637/jss.v100.i05 https://doi.org/10.18637/jss.v100.i05.
• Bürkner P, Gabry J, Kay M, Vehtari A (2024). “posterior: Tools for Working with Posterior Distributions.” R package version 1.6.0, https://mc-stan.org/posterior/. Vehtari A, Gelman A, Simpson D, Carpenter B, Bürkner P (2021). “Rank-normalization, folding, and localization: An improved Rhat for assessing convergence of MCMC (with discussion).” Bayesian Analysis.
• Chang W (2023). extrafont: Tools for Using Fonts. R package version 0.19, https://github.com/wch/extrafont.
• Eddelbuettel D (2024). digest: Create Compact Hash Digests of R Objects. R package version 0.6.37, https://CRAN.R-project.org/package=digest.
• Eddelbuettel D, Francois R, Allaire J, Ushey K, Kou Q, Russell N, Ucar I, Bates D, Chambers J (2024). Rcpp: Seamless R and C++ Integration. R package version 1.0.13-1, https://dirk.eddelbuettel.com/code/rcpp.html, https://github.com/RcppCore/Rcpp, https://www.rcpp.org. Eddelbuettel D, François R (2011). “Rcpp: Seamless R and C++ Integration.” Journal of Statistical Software, 40(8), 1-18. doi:10.18637/jss.v040.i08 https://doi.org/10.18637/jss.v040.i08. Eddelbuettel D (2013). Seamless R and C++ Integration with Rcpp. Springer, New York. doi:10.1007/978-1-4614-6868-4 https://doi.org/10.1007/978-1-4614-6868-4, ISBN 978-1-4614-6867-7. Eddelbuettel D, Balamuta J (2018). “Extending R with C++: A Brief Introduction to Rcpp.” The American Statistician, 72(1), 28-36. doi:10.1080/00031305.2017.1375990 https://doi.org/10.1080/00031305.2017.1375990.
• Gabry J, Mahr T (2024). “bayesplot: Plotting for Bayesian Models.” R package version 1.11.1, https://mc-stan.org/bayesplot/. Gabry J, Simpson D, Vehtari A, Betancourt M, Gelman A (2019). “Visualization in Bayesian workflow.” J. R. Stat. Soc. A, 182, 389-402. doi:10.1111/rssa.12378 https://doi.org/10.1111/rssa.12378.
• Grolemund G, Wickham H (2011). “Dates and Times Made Easy with lubridate.” Journal of Statistical Software, 40(3), 1-25. https://www.jstatsoft.org/v40/i03/.
• Hartig F (2024). DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models. R package version 0.4.7, http://florianhartig.github.io/DHARMa/.
• Hester J, Wickham H, Csárdi G (2024). fs: Cross-Platform File System Operations Based on ‘libuv’. R package version 1.6.5, https://github.com/r-lib/fs, https://fs.r-lib.org.
• Küng P (2023). wbCorr: Bivariate Within- and Between-Cluster Correlations. University of Zürich. R package version 0.1.22. https://github.com/Pascal-Kueng/wbCorr.
• Lüdecke D, Patil I, Ben-Shachar M, Wiernik B, Waggoner P, Makowski D (2021). “see: An R Package for Visualizing Statistical Models.” Journal of Open Source Software, 6(64), 3393. doi:10.21105/joss.03393 https://doi.org/10.21105/joss.03393.
• Müller K, Wickham H (2023). tibble: Simple Data Frames. R package version 3.2.1, https://CRAN.R-project.org/package=tibble.
• Pedersen T (2024). patchwork: The Composer of Plots. R package version 1.3.0, https://github.com/thomasp85/patchwork, https://patchwork.data-imaginist.com.
• R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.
• Schauberger P, Walker A (2024). openxlsx: Read, Write and Edit xlsx Files. R package version 4.2.7.1, https://CRAN.R-project.org/package=openxlsx.
• Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.
• Wickham H (2023). forcats: Tools for Working with Categorical Variables (Factors). R package version 1.0.0, https://CRAN.R-project.org/package=forcats.
• Wickham H (2023). stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.5.1, https://CRAN.R-project.org/package=stringr.
• Wickham H (2024). rvest: Easily Harvest (Scrape) Web Pages. R package version 1.0.4, https://CRAN.R-project.org/package=rvest.
• Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to the tidyverse.” Journal of Open Source Software, 4(43), 1686. doi:10.21105/joss.01686 https://doi.org/10.21105/joss.01686.
• Wickham H, Bryan J (2023). readxl: Read Excel Files. R package version 1.4.3, https://CRAN.R-project.org/package=readxl.
• Wickham H, François R, Henry L, Müller K, Vaughan D (2023). dplyr: A Grammar of Data Manipulation. R package version 1.1.4, https://CRAN.R-project.org/package=dplyr.
• Wickham H, Henry L (2023). purrr: Functional Programming Tools. R package version 1.0.2, https://CRAN.R-project.org/package=purrr.
• Wickham H, Hester J, Bryan J (2024). readr: Read Rectangular Text Data. R package version 2.1.5, https://CRAN.R-project.org/package=readr.
• Wickham H, Hester J, Ooms J (2023). xml2: Parse XML. R package version 1.3.6, https://CRAN.R-project.org/package=xml2.
• Wickham H, Vaughan D, Girlich M (2024). tidyr: Tidy Messy Data. R package version 1.3.1, https://CRAN.R-project.org/package=tidyr.
• Wilke C (2024). ggridges: Ridgeline Plots in ‘ggplot2’. R package version 0.5.6, https://CRAN.R-project.org/package=ggridges.
• Xie Y (2024). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.49, https://yihui.org/knitr/. Xie Y (2015). Dynamic Documents with R and knitr, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, https://yihui.org/knitr/. Xie Y (2014). “knitr: A Comprehensive Tool for Reproducible Research in R.” In Stodden V, Leisch F, Peng RD (eds.), Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595.
• Zhu H (2024). kableExtra: Construct Complex Table with ‘kable’ and Pipe Syntax. R package version 1.4.0, https://CRAN.R-project.org/package=kableExtra.